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Bound and resonance electronic states in impure graphene are studied. Short-range perturbations 
for defects and impurities of the types "local chemical potential" and "local gap" are taken into 
account. Zero gap and non-zero gap kinds of graphene are considered. A qualitative analysis of the 
electronic spectrum general features and numeric calculations are presented. A dependence of the 
resonance widths on the gap value, angular momentum values and the perturbation amplitudes is 
investigated. 

I. INTRODUCTION 

The Dirac equation is one of keystones of the relativistic field theory. However, it is an important model in the 
non-relativistic solid state theory as well. Superconductors with c?— pairing [T], the Cohen-Blount two-band model 
of narrow-gap semiconductors Q , Q , electronic spectrum of the carbon tubes form an incomplete list of the non- 
relativistic applications of this equation. During the last two years extremely much attention was payed to the problem 
of the electronic spectrum of graphene (see a review ^ ) . Two-dimensional structure of it and a presence of the cone 
points in the electronic spectrum make actual a comprehensive study of the external fields effect on the spectrum and 
other characteristics of the electronic states described by the Dirac equation in the 2-1-1 space-time. We consider in 
this work the resonance and bound states of the 2-1-1 Dirac equation due to the short-range perturbation using some 
approaches developed in Ref.[5,] .We do not take into account the inter- valley transitions. Particular attention to this 
case stems from the effectiveness of short-range scatterers in contrast to the long-range ones: an effect of the latter is 
suppressed by the Klein paradox 0. Our work takes into account the obvious fact that the Kohn-Luttinger matrix 
elements of the short-range perturbation calculated on the upper and lower band wave functions are not equal in a 
general case. This means that not only the potential but the mass perturbation can be present in the perturbed Dirac 
equation. Resonance states in the zero-mass Dirac model of graphene were considered for distinct perturbations in 

a, i, a. 

The paper is organized as follows. In the first section we formulate the problem in terms of the 2-1-1 Dirac equation 
with the delta function perturbation taking into account a local change of the potential and mass (gap) due to the 
crystal defect. An exact solution of this problem leads us to the characteristic equation determining the bound and 
resonance states in graphene. In the second section we analyze the characteristic equation for the zero-gap graphene 
case. Both local and global behavior of resonance states for various magnitudes of the perturbations are investigated. 
Analytic and numeric results are presented. Just analytic approach allowed us to investigate the case of exponentially 
narrow resonance that could be difficult to do numerically. In the third section we consider the bound and resonance 
states in graphene with non-zero mass (gap). The exponential approach of the bound state energy levels to the band 
edge was studied analytically; numerical analysis is presented as well. 

II. CHARACTERISTIC EQUATION 

The Dirac equation describing electronic states in graphene reads [3] 



(1) 



where vp is the Fermi velocity of the band electrons, 7^ are the Dirac matrices 

7o = c^3, 7i = c^i, 72 = 
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CFi are the Pauli matrices, 2touf^ = Eg is the electronic bandgap, ip (r) is the two-component spinor. The electronic 
gap can appear in the graphene monatomic film lying on the substrate because of the sublattices mutual shift p^ : 
the spin-orbit interaction can be the reason too. The spinor structure takes into account the two-sublattice structure 
of graphene. 5m (r) and ^(r) are the local perturbations of the mass (gap) and the chemical potential. A local mass 
perturbation can be induced by defects in a graphene film or in the substrate jlOl] . We consider here the delta function 
model of the perturbation: 

(5m (r) = — 65(r — ro), V{r) = —a5{r ~ ro), (2) 

where r and tq are respectively the polar coordinate radius and the perturbation radius. Such short-range perturbation 
was used in the (3+l)-Dirac problem for narrow-gap and zero-gap semiconductors in 0]. 
The perturbation matrix elements 

diag{Vi,V2)5{r ~ r^) (3) 

are related to the a, b parameters as follows 

Vi=-{a + b), V2 = b-a (4) 

The delta function perturbation is the simplest solvable short-range model. Finite radius tq plays a role of the 
regulator and is necessary in order to exclude deep states of the atomic energy scale. The finite perturbation radius 
To leads to the quasi-momentum space form-factor proportional to the Bessel function that justifies our neglect 
of transitions between the points K and K' Ref. _3]. The two-dimensional Dirac problem with the scalar short- 
range perturbation Eq. ([2|) (but without the mass perturbation) was considered in Ref. [13]. The obtained there 
characteristic equation for the discrete energy spectrum contains one mistake. We corrected it in our previous work 
Ref. 12] and took into account the mass perturbation 6m (r) . Here we present more detailed analysis of the electronic 
bound states and new results on resonance states in the zero-gap and non-zero-gap graphene. Particular attention is 
paid to the effect of relative intensities of two perturbations Vi and V2 (a and b). 

Let us present the two-component spinor in the form 



exp {—iEt) 



gj {r) exp [i {3 + 1/2) Lp] 



where j is the pseudospin quantum number; j = ±1/2, ±3/2, .... In opposite to the relativistic theory, this quantum 
number has nothing to do with the real spin and indicates a degeneracy in the biconic Dirac point. The upper fj (r) 
and lower Qj (r) components of the spinor satisfy the equations set 

^ + -^g, -{E-m) f, = {a + b) 5{t - ro)f„ (6) 
- 1^ + -iE + m) g, ^{a-b) d{r - r,)g,. (7) 



These equations have a symmetry: 



fo ^ 5j> E -E, j ~j, a -a. (8) 
Let us introduce the function ipj (r) = fj/gj- It satisfies the equation: 

^-%,^E{ip'^+l) /[ia + b)^] + ia-b)]+S{r-ro)=0 (9) 
Integrating in the vicinity of r = tq 

e^oj^^^ro-s) (a + 6)(p2 + (a-6) 

we obtain the matching condition 
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arctan (^ip^ ^{a + b) / {a - b)j - arctan (^ip^ ^(a + 6) / (a - b)^ = \/ a? -b"^^ c? > 6^, (11) 

. The upper and lower component matching conditions resuhing 



where ip^ = p^ (ro - 6) , ipj = p^ {vq + 5) , 6 
from Eq. (fTT|) read 



where 



A 

A 



cos y/a^ — b'^, 

a+b 



'^sinVa2 -62 



f^sinV^^^, 



a'^b'>0 



cos Va^ — 6^ 

is the orthogonal rotation matrix. It transmutes into the orthogonal boost matrix for 6^ > 



(12) 



(13) 



A 



cosh \/b'^ — a^, 



'|±2- sinh - a2 



sinh \/W~~a^, cosh 



> 0. 



(14) 



The general solution can be found solving the second-order equation obtained by excluding one of the spinor compo- 
nents from the equation set Eq. Eq. ([7]) in the domains < r < rg and r > tq : 



III 



'j^2 2 j U ~ 1) 
h — m 



fj = 0. 



(15) 



This equation is related to the Bessel one. We assume at first E to be real and satisfying the inequality < m^. 
Then the general solution of Eq. (fT5)) reads 



fj = CiV^/j_i/2 (Kr) + C2VrKj_i/2 (nr) , 



(16) 



where k = V m'^ — E"^ is the principal value of the root; (z) and Ki, (z) are the modified Bessel functions. The 
constant C2 = in the domain < r < tq, while Ci = in the domain r > tq. It is useful to introduce the following 
notations: 



Kj (z) = ifj_l/2 {z) /Kj+i/2 {z) , Xj (z) = /j_i/2 [z) /Ij+i/2 {z) 

Expressing the gj-component from Eq. ([7]), we can write 



= ^{m + E) / {m- E)Ij {nro) , 



(17) 



(18) 



p+ = ^{7n + E)/{m- E)JCj (^m) 



(19) 



Substituting the expressions Eq. p^ . Eg. p^ into the matching condition Eq. (dJ), we obtain the characteristic 
equation for the bound state energy levels (for E"^ — < 0): 



K \Ij (nro) — ICj {nro)] — T (a, b) [(m — E) {a ~ b) + {a + b) (m + E)!^ (k^q) JC^ (k^q)] 
where T (a, b) .is determined as follows: 



(20) 



T{a,b) = 



tan {^/a^^ - b^) /Va^ - b^ if > 6^ 



tanh (Vfe2 - a2) if b^>a^. 

We presented above an analysis for the case of > b^; the case of b^ > c? can be considered similarly. 



(21) 
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This equation turns to the characteristic one obtained in for 6 = apart from the mistakenly omitted terms 
in the right hand side of Eq. (PU)) . This characteristic equation is unambiguously determined for the bound states 
with energy levels lying in the real axis segment [— m, m] . 

We write this equation in another form making the symmetry Eq. ([H]) manifest: 



= T(a,6) [(m - £:) (a - 6)/j+i/2 (Kro) (k^o) + (a + 6) (m + S)/j_i/2 (Kro) Xj_i/2 (Kro)] , (22) 

This equation was derived for a study of bound states situated in the gap that will be done in section IV. 

An analytic continuation of this equation from the energy real axis segment {—m, m) onto the energy complex 
plane is necessary in order to study resonance states. Using the obvious relations 

-ip, if i? > TO 



k = \/ — 



ip, if E < — m, 



where p — \AE^^^TO?,and the known relations between the Bessel functions [K 



exp {—i-Kv 12) {z exp (i7r/2)) , — tt < arg z < 7r/2, 

exp (z7ri^/2) (2; exp (— i7r/2)) , — 7r/2 < argz < tt, 

Z7r/2exp {iiru /) Hi^^ (zexp (i7r/2)) , — tt < argz < 7r/2, 

r(2) 



—177/2 exp {—iTTv/2) Hi [z exp {—111/2)) , — 7r/2 < argz < tt, 
we obtain the characteristic equation in the form 



P 



Jj (pro) + (pro) = -T {a, b) (E - m) {a - b) - {E + to) (a + b) J, (pro) H)"' (pro) 



(a) 



(23) 
(24) 

(25) 



where 



a 



1 if Re > TO, 
2 ifRci;<-m. 



We have introduced the notation: 



J, (Z) = J,-l/2 (Z) / J, + l/2 (Z) , (Z) = ijj!) /2 (Z) /ijjj,/, (Z) . 



(26) 



The cases of = b^ (coordinate angles bisectrices in the (a, 5)-plane) are degenerate: the energy imaginary part 
vanishes in this limit since according to Eq. (|4]), one of the matrix elements Vi(2) equals zero. For instance, the 
effective second-order equation for the upper spinor component takes the simple form: 



E^ -m^ 



J (j - 1) 



+ {E + m) Vi (r) 



= 0. 



(27) 



This equation obviously does not contain a resonance state because of absence of terms of the type characteristic 
for the relativistic resonances V^, dV/dr etc [14 j. 



III. ZERO-GAP GRAPHENE 



Let us consider Eq. ((25)) in the zero-gap case to = 0. The symmetry Eq. ([8]) allows us to restrict the analysis by 
the right energy half-plane. If we consider explicitly only the right energy half-plane (a — 1), we should take into 
account both positive and negative angular momentum quantum numbers j = ±1/2, ±3/2, . . .and both positive and 
negative potential amplitude a values in order to obtain a complete picture, while considering both a = 1 and a = 2 
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cases it was enough to take into account only positive j. The characteristic equation can be essentially simplified in 
the zero-gap limit m = : 

E { J,_i/2 {Ero) Hf^^^^ (Ero) + J,+i/2 (Ero) Hf\^^ (Ero) + 
- T (a, b) [{a + b) J,_i/2 (Ero) Hf\^^ (Ero) - (a - 6) J,+i/2 (Ero) H^%/^ (Ero)] } = 0, (28) 

j = ±l/2,±3/2,- • • , ReE>0. 

Notice that this equation has a root E ^ 0. Vanishing of the energy imaginary part stems obviously from vanishing 
of the free Dirac density of states g {E) — TrlmGr {E) at E ~ 0. Here Gr {E) is the free retard Green function. 

Complex roots of the characteristic equation will be interpreted here as resonance states. Now we will study a 
distribution of the complex roots of Eq. P5|) . i. e. resonances. 

It is convenient to divide the plane (a, b) into four regions separated by the coordinate angles bisectrices a? = b^ 
(see Fig. 1). We numerate these regions from I to IV. 

A. Resonances: a? > 

Let us consider at first the region I of the (a, 6)-plane (see Fig. 1). It is useful to introduce there the hyperbolic 
variables: 

a = Picoshtpi, 6 = p]^ sinhf/)]^. (29) 

Then Eq. ^ takes the form: 

J,_i/2 (Ero) i?j;\/2 (Ero) + J,+i/2 {Ero) H^,% (Ero) 



tan pi 



exp (V'l) Jj_i/2 (Ero) Hf\^^ (Ero) - exp (-^^i) J,+i/2 [Er^) H^%/^ (Sro) I . (30) 



Periodicity of tanpj^ allows us to restrict the analysis here by the segment < pi < tt. Consideration of only right 
half-plane of energies Re i? > and j = ±1/2, ±3/2, . . .will be suffice due to the symmetry Eq. 

For a fixed pi, we can consider the asymptotic behavior at ~^ oo. Eq. (j30p asymptotically transforms into the 
equations 

Jj_i/2 (Ero) Hf\^^ (Ero) = 0, Vi ^ oo, (31) 



Jj+i/2 (Ero) i/j^\/2 {Ero) - 0, ^ -oo. (32) 

Since the Hankel functions ^^j^\/2 (^''o) have no roots in the right halfplane [l^, solutions £',V^^''^-'of Eq. pTjl . Eq. 
([5^ read 

^^i'"'^'^ =A(/"-i/2)/ro, Vi^oo, n=l,2,... (33) 



i?(/+i/2) = A^f+i/^V'^o, V'l^-oo, n=l,2,... (34) 

with Xl^^^^'^^ being the roots of the Bessel functions Jjipi/2 (z) , which are real and tend to infinity. Thus we conclude 

that solutions of Eq. ([50]) are distributed on the arcuated curves beared on the real axis in the points En^^^"^^ (see 
Fig. 2) 
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Notice that the case of oo corresponds to ^ (see section II, Eq. (|27l) . When we move along the 

non-singular hyperbolae € (0,7r/2) from V'l — ~oo to = +00, the n-th arc-like curves are circumscribed in 

the complex plane of E from the lying on the real axis point En^^^^'' = rg^^A^^^^^^'to the similar point En = 
Xl^^^^^\ When approaches n/2 (singular hyperbola), the "arc" height tends to infinity. Moving along the 

arc-like curve takes place anti-clock-wise, when ipi varies from —00 to +00. A transition between the points 
corresponds to swapping of the upper and lower components of the wave function spinor amplitudes. However, crossing 
the singular hyperbolae 7r/2 leads to an essential rearrangement of the arc structure. Notably, the function tanpj^ 
changes its sign and alternative combinations of nearest neighbor Bessel's function roots are connected by arcs in 
order to satisfy the characteristic equation. Thus, the resonance arc-like curve begins in the point and moves 

clock-wise into the point E^^^^^'^\ when 7r/2 < p^ < it . K schematic picture of such transitions of arcs is presented 
in Fig. 3. Results of numerical calculations for p^ g (0,7r) presented in Fig. 4 for pi S (0,7r/2) and in Fig. 5 for 
Pi £ (7r/2, tt) . They are in good qualitative agreement with our analytic consideration. A periodicity with the period 
lS.pi = TT takes place. 

Motion along rays in the region I of the (a, 6)-plane maps into closed resonance trajectories in the energy complex 
plane in the case of the ray lying in the upper half-plane. In the case of the ray lying in the lower half-plane we 
obtain repeating jumps from one resonance to the next one instead of the closed curves (see Fig. 6).. When the ray 
approaches the coordinate angle bisectrix, the closed curves diameter decreases so that the curves shrink into a point. 

Consider now the left quarter (region III in Fig. 1): 

a = — P3 cosh-^a, b — —p^smhij}^. (35) 

Then the characteristic equation reads 

J,_i/2 {roE) Hf^^,^ {roE) + J,+,/, {roE) H^'\^^ {r,E) 

= tanp3 cxp(-V'3)'^j+i/2(^o-E)i?j^J^\/2(''o-B)-exp(V'3) Jj_i/2(ro£;)i?]i\/2(7-o-B) . (36) 
The asymptotic form of the equation at jV'al — >■ 00 is following 

Jj-i/2 (roE) = 0, V3 00, 



Jj+1/2 (roE) = 0, tp3^ -00. 

When 7/13 varies from —00 to +00 circumscribing the hyperbola in the region III, we have motion from the roots 
of Jj_|_i/2 (foE) to the roots of Jj-1/2 {tqE) in the energy complex plane (the same arc-form curves circumscribed 
anti-clock-wise or clock-wise respectively for pg G (0, 7r/2) or (7r/2,7r)). Numeric calculations are similar to the case 
of the region I; they confirm conclusions of our analysis. 

An increase of the potential radius rg shifts the roots in the direction of the coordinate origin. 

B. Resonances: > 

Now we consider the region II in the (a, 6)-plane (see Fig. 1). We introduce the hyperbolic variables here differently: 

b — p2COship2, a = p2Smhip2- (37) 

The characteristic equation takes the form 

J,-i/2 (roE) H^%/^ {roE) + J,+i/2 (roE) Hf\f^ (r^E) 

exp (V'z) Jj-1/2 (roE) Hf\/^ (roE) + exp (-1^3) Jj+1/2 (roE) i?]^\/2 (ro^)] tanhp. (38) 
One can see that in the limit of a ^ b we have the equation 

Jj-1/2 (roE) =0, -02 ^ 00, 
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while for a ^ —b we have 

-^5+1/2 (roE) =0, tp2^ -oo. 

Thus, the roots chcumscribe arc- form curves cfock-wise from the point A^^^^^-'/ro to A^^^^^^V^o (see Fig. 7), while 
-02 varies from +oo to —oo (motion along the upper hyperbolae in Fig. 1 anti-clock- wise). 

Let us learn a limiting form of the arcs, when — * co- As tanh p2 — *■ 1, Eq. ()38|) takes the form 

[J, (roEjn) - exp (-^3)] [Hj [roEj^] - exp (-^-2)] = 0, (39) 

where Ejn = Ejl^ — iV jn/2 is the complex root of the characteristic equation; we used the standard notions for the 
real and imaginary parts of the resonance energy. 

Notice that roots of the second factor are not physical and must be ignored. Eq. ([55]) can be written in the form 

Jj (^0 (Ef^ - «r,„/2)) = exp , P2 00, (40) 

Since Jj {z) is real on the real axis, we have Jj {z)* = Jj {z*) and, therefore, Jj {^0 {^E^^^ ~ iV jn/^j^ = 

Jj {tq [sf^ -|-irj„/2^^ = exp(-V'2) ■ So we see that E^^} + iTjn/2 is a root of Eq. (gD]) as well. Since the 

ends of the arc are simple zeros of the Bessel function, these two arcs coincide and are flat in this limit. Eq. (HD|) 
takes the form 



J, (ro£;j,V) = J,-i/2/ J,+i/2 = exp (-V2) = ^/{h-a)l{h + a), 6^ » ^2 _ (4^) 

We see that the energy imaginary part is getting extremely small and the arcs approach the real axis segments 
[Afj-'~^^^V''o, A,^^^^^^V^o] so fast at P2 — > 00 that it is difficult to observe a presence of the imaginary part numerically. 
However, this process can be followed analytically. In order to investigate the limiting behavior of the energy imaginary 
part, we expand Eq. (I38p for p2 — > 00 taking account of the asymptotics leading terms. In particular, the expansion 
tanhx ^ 1 — exp (—2a;) -|- . . . for a; ^ 1 was used. Taking the real part of Eq. ([38|) . we obtain after some algebra 

ReJjiroE,,,) [1 - (1 - exp(-2p2) exp (7/^2) ^eUj {tqEj^))] 
= - Ren, {roEjn) + TnJ- {roE^n^ (ro/2) (1 - exp (-2^2)) exp (-^2) Im (roS,„) 
+ (1 - exp (-2P2)) exp (--02) . 
Expanding this formula for small Tjn we obtain in the leading order 



exp(-V.2-2fe) ReH,(roi;W)+exp(-V.2) 
J, [r,E^^ Imn, [roEll^) 

Here Ejl^ is determined from Eq. (j4ip . It is seen from Eq. (|42|) that the resonance width Tj„ decreases exponentially 
at large pj- 

Returning to the initial notions we find: 



roTjn = 2exp(-2v/62 



V{b-a)/{b + a)Jf+if2 



Re 


"rr(l) /rr(l) 
-"j-l/2/-"j + l/2 




a) / {b + a) 


Im 


"rr(l) /ttW 

^j-l/2/^j+l/2 





b' » a' 



(43) 



The obtained formula indicates that these resonances are extremely sharp for the case of small \j\ and, therefore, 
they can give a strong resonance scattering of electrons. Using the Bessel functions J^, asymptotics for large v 
[ist , we obtain the following estimate for the resonance width asymptotic at large | j | : 



(r-i) v/(6-a)/(6 + a)exp [-2^b^ - a^) (|j| /e) 



(44) 
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This function catastrophically increases beginning from \j\ ~ 4. It means that the imaginary part of the energy is 
large for any reahstic P2 so that solutions with j > i are not really sharp resonances and they will not contribute into 
the electron resonance scattering. 

Results of numeric calculation for resonances in the case of > are presented in Figs. 7 and 8. Approaching 
of arc-like resonance curves to the real axis segments, when p2 increases, is shown in Fig. 7. Fig. 8 shows the 
resonance trajectory corresonding to a motion along the ray in the (a, 5)-plane region II. When the radius increases, 
the trajectory approaches a point on the real axis. 

The region IV (lower quarter) can be analyzed similarly using the hyperbolic variables: 

b — —p^^coship^, a — —p^sinhil;^ (45) 

An analisis of this region could be done similarly to the region II, but we omit it here since no qualitative distinction 
is observed in this case. Numerical results for the region IV are omitted by us as well because of this similarity. 

So we have discussed a local behavior of the resonance trajectories within the chosen quarter in the (a, 6)-plane. 
Now we will consider the global behavior, when all quarters are passed along hyperbolae with the fixed hyperbolic 
radius p. This process creates pairs of arc-form curves bearing on two common basing points lying on the real axis. 
These figures of eight are arranged at the E complex plane symmetrically relative to the imaginary axis in the case 
of m = 0. Resonances move anti-clock- wise from the root A^/^^^^-'/ro to the root A^-'~^^^Vro circumscribing the arc, 
when (j) £ [— 7r/4, 7r/4] When (p £ [7r/4, 37r/4] , a different arc appears. It connects the same points A^^'*'"'^/^-' and 
^n~^^^\ but moving in the opposite direction (clock- wise) . We obtain in result the figures-of-eight in the £'-complex 
plane. This picture is repeated with the period A0 = tt. 

Such excursion around all quarters of the (a, 6)-plane can be carried out circumscribing the circles in the (a, 5)-plane 

a = pcos(j), b^psincj). (46) 

One quarter after another is passed in this case too, but now the figures-of-eight do not touch the real energy axis 
anymore (see Fig. 9), i. e. all resonances belonging to this set have rather large widths. Numeric analysis confirms our 
prediction on the mapping of the circular excursion in the (a, b) plane onto the energy plane. When the perturbation 
intensity p is increased, the picture is getting more complicated: all the figures-of-eight are associating into a single 
whole (see Fig. 9 for large p). Notice that the resonance trajectory self-crossing points are fixed points relative to p 

changes. !They correspond to resonance set, which do not change at motion along the ray a — b. 

IV. GAPPED GRAPHENE 

While the pristine graphene is gapless, a contact with the substrate can induce some narrow gap [T^. That is why 
we study here bound and resonance states in the gapped graphene as well. 

Let us consider the characteristic equation for the case oi m ^ 0. We introduce the dimensionless parameters: 

e = E/m, rom = ro/rc = 7, (47) 

where Tc = 1/m is the "Compton radius" {h = 1, vp — !)• The characteristic equation for bound states Eq. ((22l) 
takes the form: 

Vl-e^ [lj-l/2 (z) Kj + i/2 {z) - Kj_i/2 [z) Ij + i/2 {z)] 

= T(a, b) [(1 - 6) (a - b) J,+i/2 (z) if,+i/2 [z) + (1 + e) (a + b) I^_,/2 (z) K^_^/2 {z)] , (48) 

where z = 7\/l^-?. The perturbation radius rg value is of order of the lattice spacing d, while the " Compton" radius 
is rather large: >> d : h/ (mQVp) ^ 10^* cm even for the bare electron mass ttiq; the real effective mass at the 

band extrema can be estimated as 10~^ — 10~^ of uiq. Therefore, the actual value of 7 is less than unity. 

We consider here some general properties of the bound states electronic spectrum resulting from the characteristic 

equation (|48|) and the effect of non-vanishing gap on the resonance states. Results of the numerical solution will be 

presented too. 

A. Bound states: > 

Let us assume at first j = 1/2. We consider Eq. for the energy situated in the gap < 1. We use the Bessel 
functions limiting forms for the states lying near the gap edges 13] 
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lAz)^{z/2r Xo(z)~-logz, KAz)-\t{v){z/2)- 



(49) 



in order to transform the characteristic equation ([^5]) . Here F {v) is the gamma function. We obtain a simple relation 
for small z. 



1 + (7V4) (1 - e") log (7/ 



7- 



tan \J a} — \P- ' 



1 -e2 

(1 - e) (a - b) /4 - (1 + e) (a + b) log (7/ 



1 



(50) 



Formulae like Eq. (jSOp determine a map from the plane (a, 6) onto the complex energy plane; in the case of bound 
states the map is carried out onto the real energy axis segment — 1 < e < LNow we consider the asymptotic energy 
level behavior near the upper edge of the gap e — > 1. Eq. (|50p takes the asymptotic form: 



tan y/a? 



62 



(a + 6) log [272(1-6)] 



(51) 



Therefore, we can write 



exp 



(1/7) 



1 



V^2— 52 



a + b tan Va^ - 6^ 



(52) 



It is clear that e — > 1 if a+& — > +.0. This result conforms the well known general property of the two-dimensional 
quantum systems: a threshold for creation of the bound state is absent for j = 1/2. 

Formula ([52|l can be rewritten in the hyperbolic variables in the region I (see Eq. ([29| as follows 



1 - e = (27^) ^exp[-l/ (7exp(V'i)tanpi) 



nn + (where n = 1, 2 . . .), i. e. the e 



(53) 
1 "fronts" 



It is seen from the last formula that e — > 1, when pi 
approach the hyperbolae pi — \/ a? — b"^ — *■ 7rn + 0. 

Let us study now a behavior of the bound state trajectory near the lower band edge. Using Eq. (f50|) we conclude 
that when e — > —1, the following relation holds: 



l?^7tan(pi)exp(-Vi)(l-e)/4. (54) 

This equation shows, how the bound states trajectory approaches the lower edge of the gap. Putting e = — 1 we 
obtain the equation, which determines the e = — 1 boundary lines in the (a, 6)-plane, separating the domains of the 
bound states: 

tanp]^ ~ 2expV'i/7. (55) 

When -01 — ^ +00 along one of these lines, p^ — > 'k 12 + 7rn — 0, rt = 0, 1, . . .. The line n = crosses the abscissa axis 
in the point (see Eq. ([M)) ): 

&o = 0, ao = arctan (2/7) . (56) 

Then this curve crosses the bisectrix a + 6 = in the point ai = I/7, b\ = — 1/7. Any point situated in the domain 
restricted at the left by the bisectrix a — = and at the right by the line Eq. ([55)1 (e = — 1 "front") (which 
asymptotically approaches the hyperbola p^ = tt/2 — 0, tpi — > +00), gives a solution of the characteristic equation 
lying in the gap |e| < 1 that is an eigenvalue of our problem. A motion from the bisectrix a + & = to the curve Eq. 
(1551) in the region I of the (a, 5)-plane is mapped to the motion from the upper to lower edges in the energy gap. 

Other domains with n = 1, 2, ... , lie between approaching the hyperbolae pi = im — curves (e = 1 "front", 
Itpil — > 00) and the e — —1 "front", approaching the hyperbolae pi = tt/2 + im — at -01 — +00 and the 
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hyperbolae pi = Trn + at -01 — > — oo. These Unes boarder the countable set of domains in the (a, 6)-plane, where 
the bound states exist. When we move along the lying in the region I ray beginning in the coordinate origin in the 



direction of the pf 



b increase, we cross this set of domains so that every time the bound states trajectory 



starts at e = 1 and terminates at e = — 1. 

A dependence of the bound state energy on the radius pi, at the motion along the rays in the (a, 5)-plane (see inset 
to Fig. 10) was investigated numerically in the case of > 6^. The motion between the coordinate origin and the 
first e = — 1 "front" maps into the mildly sloping transition between the e = 1 and e = — 1 edges of the gap (solid line) 
and the sharp sloping transition respectively for the rays in the upper and lower half-planes of the (a, 6)-plane (first 
pair of the curves). The second and so on pairs of curves show similar transitions between the domain boundaries. 

Now we consider the region III. Let us return to Eq. ([50]) that is valid near the band edges and everywhere inside 
the gap if 7 << 1. Using the hyperbolic variables Eq. (|55|) we can write for the bound states near the upper edge 
e = +1 : 



1 - e - (2Y 



exp 



exp (V'g) 



7tanp3 



(57) 



Eq. (|Ff|) is satisfied near the hyperbolae p^ = irn — (then tanpg — > —0). This determines the e = 1 "front" of 
domains containing the eigenvalues. The e = I "fronts" are formed by the approaching the hyperbolae p^ — nn ~ 0, 
(n = 1,2,.. .) at IV'al — > +00. Therefore, we have the exponential approaching of the energy level to the upper band 
edge in the region III similarly to the region I. 
The e = -1 "fronts" Eq. jMl) 



approach the hyperbolae 
n = 1,2,.... 



nn 



tanp3 = - (2/7) expo's 
— 0, when ip^ — > —00, and hyperbolae p^ 



-7r/2 + nn + 0, when ip^ > +00, 



B. Bound states: b > a 

Let us consider the region II in the (a, 6)-plane.The approximate characteristic equation Eq. ([50)1 for the states 
near the band edges in the case of 6^ > reads: 



1 + (7V4) (1 - e2) log (7/ 



l-e- 



7- 



tanh y/b'^ — a'^ 



y/b^ - a2 

We obtain for the energy near the upper band gap edge 



(1 - e) (a - b) /4 - (1 + e) (a + b) log (7V1 - e^) 



(58) 



l-e'^ (27^) ^ exp 



(I/7) 



1 



a + 6 tanh \/b'^ - 

Thus, the e = 1 "front" approaches the bisectrix a + 6 — > exponentially. We introduce here the variables 

6 = /92 cosh ?/'2 , a = sinh ?/;2 ■ (59) 

Then Eq. ([58l) takes the form 



1 + (7V4) (1 - e') log (7/ 



-7tanhp2 (1 - e) exp {-tp^) /4 + (1 + e) exp [ip^) log (^7 V 1 - e^) 
We obtain for the energies near the upper band edge e — > 1: 



(60) 



1-e^ 1/ (272)exp 



exp (-i/'a) 



7 tanh P2 



(61) 
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Thus, the e = 1 "front" approaches the bisectrix a + 6 = at — * 0. There is no "front" at the bisectrix a — 6 = 
as it was shown above. Eq. (|60p can be rewritten for — * oo ; 



[1 + (7/4) (1 - e) cxp i-ij^)] [1 + 7 (1 + e) exp (V-a) log (7^1^) 



= 0. 



Therefore, 



7(l + e) exp('02)log (-fV 



1 



-1. 



(62) 



(63) 



Eq. ((63|) is satisfied if e — > —1 and -02 — * +00. Therefore, the e — —1 "front" approaches the bisectrix a — b ~ 0, 
if P2 — > 00. This "front" hes in the region I (see above). Thus, a motion in the region II of the plane (a, b) from the 
bisectrix a + 6 = to the bisectrix a — b = corresponds to a motion of the energy eigenvalue from the upper band 
edge e = 1 to the lower band edge e = — 1, when p2 >> 1. 

A dependence of the bound state energy on the radius P2, at the motion along the rays in the (a, 5)-plane (see inset 
to Fig. 11) was investigated numerically in the case of < 6^. The motion along the ray maps into the monotone 
transition of the energy level from the band gap edge e = 1 to the constant situated in the gap. Motion along the ray 
in the right half-plane maps into the faster energy dependence and with lower asymptotic energy value, than in the 
case of the ray in the left half-plane. When the ray approaches the bisectrix a — 6 = 0, the asymptotic value tends to 
e = -l. 

Let us consider now solutions of the characteristic equation near the band edges for the parameters a, 5 in the region 
IV. We introduce the hyperbolic variables (see Eq. ([45])). 
Then Eq. ([55)1 can be approximately written in the form: 



7tanhp4 (1 - e) exp (-^4) /4 + (1 + e) exp (-l-V'4) log (^^\/ 
1 = (7/2)exp(-V'4)tanhp4 



The line 



1-e 



(64) 



(65) 



corresponds to the e = — 1 "front" in the region IV. This line is a continuation of the line Eq. (|55p : they meet in the 
point a — —b — I/7. The line Eq. (j65|) turns into the asymptotic at — > 00 : 



■04 

Let us transform Eq. in the limit of p^ — > 
[1- (7/4) (1-6) exp (-^4)] 
The "front" e — > 1 equation takes the form 

1 - 



-log (2/7). 



7(l + e)exp (-04) log (7^ 



= 0. 



(66) 



(67) 



(4/7) exp (V'4) 



(68) 

Thus, the "front" e — *■ 1 approaches the bisectrix b + a = asymptotically, when ■04 — > —00. Therefore, a motion 
from the line Eq. (j66|) to the bisectrix & + a = in the (a, b) plane is mapped to a motion from the lower edge of the 
gap to the upper one when p^ » 1. 

If 7 increases, we conclude from Eq. ([55)1 that the second front 1 -I- e — > moves to the bisectrix 6 — a = 0, i. e. 
an increase of the mass widens the domain of the region IV, where the eigenvalues exist; it fills all region IV in the 
extreme limit. of large 7. Notice that this case is not realistic for an impurity problem, but can be important for the 
quantum dot case Q. 

When 7 << 1, it is seen from Eq. that the "front" l + e — > tends to the "front" 1 — e — > at V'4 — > —00, 
i. e. the fronts tend to merge, when the gap tends to zero. 

Thus, our analysis shows that there is a countable set of eigenvalue domains in the case > b^. The boundaries 
of these domains approach the hyperbolae described above. In the case of b^ > a^, we have a saturation instead of 
periodicity, when radius increases. 
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C. Bound states for higher angular momenta 

Let us consider now the eigenvalue spectrum for the angular momentum j = 3/2, for instance, when the parameters 
a, b lie in the regions I and II.. Using the Bessel functions expansion for small arguments Eq. (|48p can be written 
as follows: 

1/7 - 7 (1 - e') /12 = T {a, b) [{a - fe) (1 - e) /12 + {a + b) {1 + e) /4] . (69) 
In the region I, this equation takes the form 

1/7 - 7 (1 - e') /12 - tanpi [(1 - e) exp {-ib,) /12 + (1 + e) exp (^1) /4] . (70) 
The "fronts" e — > 1 have the form 

tanpi = (2/7)exp(-Vi)- (71) 



We have asymptotically: 



tt/2 + vrn — for ipi — > —00, 
pA^i)^{ n-1,2,... (72) 

7rn + for ipi — > 00, 



The "fronts" e — > —1 have the form 



Thus we have the lines 



tanpi = (6/7) exp (-^i) • (73) 



TTn + for V^i > —00, 

Pi(V'i)^<! n = l,2,... (74) 

7r/2 + TTfi — for > 00. 

The "fronts" Eq. ((7T|) (e — > 1), Eq. ([73]) (e — > —1) cross for the first time the abscissa axis respectively in the 
points 

ai — arctan(2/7), 61 = (75) 
02 = arctan(6/7). (76) 

Therefore, a motion from the point ai to the point 02 maps into a motion from e = 1 to e = —1. More than that, having 
in mind a continuous dependence of solutions on the parameters a and b the "fronts" Eq. (|7ip . Eq. (j73[) obtained 
from the asymptotic equation (I70p cross one another. This prediction is confirmed by the numerical calculations (see 
Fig. 13). Exclusively fast dependence of the energy eigenvalues on the parameters a and b will be observed in the 
vicinity of the crossing point. 

Let us consider now the region II. We have here from Eq. ((69|) : 

1/7 - 7 (1 - e') /12 = tanhp2 [- (1 - e) exp (-^1) /12 + (1 + e) exp (^1) /4] . (77) 
The "front" e — > 1 is given by the equation 

tanhp2 - (2/7)cxp(-^2)- (78) 

Since < tanh /52 — 1 iP2 ^ 0) i we have 

V'2>-log(7/2). (79) 

The domain, where the energy eigenvalues exist in the case of 7 = 2, is bounded by the positive half- axis b and 
the first quadrant bisectrix. When 7 < 2 and decreases, the "front" e = 1 turns about the coordinate origin from the 
ordinate axis approaching the first quadrant bisectrix. The eigenvalues domain is narrowing in result. When 7 > 2 
and increases, this boundary turns anti-clock- wise approaching the second quadrant bisectrix. When 'ip2 — > 00, Eq. 
(iTTl) can be satisfied only if 

H-e = 4/(7tanhp2)exp(--02): '^2 — ' °o (80) 
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We have for — ^ oo (higher hyperbolae): 

1 + e = (4/7) exp (-V-a) , i'2 (81) 

When V'2 is fixed, the asymptotic magnitude of the eigenvalue £as at p2 — > 00 can be found from the quadratic 
equation 

eL - 7"' [6 exp (^-3) + 2 exp (-^3)] 6,, - 1 + 12/7^ + (2/7) exp (-^-2) - (6/7) exp (^2) = 0. (82) 

The domain of eigenvalues in the region IV is narrower in the case of j — 3/2, than for j — 1/2 similarly to the 
region II. Generically, an increase of j at fixed 7 leads to narrowing of the eigenvalues existence domains. 

The case of higher angular momenta j = 5/2, . . . can be analyzed similarly. The energy discrete eigenvalues exist 
for all J, but in opposite to the case of j = 1/2, there is a threshold for the bound state appearing at the upper 
edge of the gap. Existence of bound states for all higher angular momenta is a specific feature of the delta function 
perturbation. 

The states with j < can be easily obtained with a use of the symmetry transformation Eq. ^ . 



D. Diagram of the bound states 

General results of mathematical and numerical analysis for a distribution of bound states for j = 1/2 and j = 3/2 
are presented respectively in Figs. 12, 13. Shaded domains represent basines of the bound states existence. The thick 
lines correspond to the upper edge of the band gap e = +1, while the thin ones represent the lower edge e — —I. 

In the case of j = 1/2, domain boundaries approach asymptotically the hyperbolae. One crossing point of boundaries 
can be seen. Therefore, we have one twisted domain in this case. Notice that this crossing looks so in a projection 
onto the (a, 6)-plane. A three-dimensional picture of this feature is shown. in Fig. 13: it is seen that the "fronts" 
e = ±1 lie at different "floors" so that the crossing curves in the plane curves appear to be skew ones. 

Notice that our diagrams of states are obtained assuming single- valued dependence of the bound state energy on the 
parameters a and b. This is true everywhere in the (a, 6)-plane except very small islands near the domain boundaries 
(see Fig. 10 of this paper and our work [l2|, where the twisting point neighborhood was analyzed numerically). 

In the case of j = 3/2, the bound state existence domains are getting narrower, a number of the boundaries crossing 
points increases. The diagram of states for higher j is given us in Ref. [l^ . 



E. Resonances: m > 



It is necessary to return now to Eq. ([251) written in terms of non-modified Bessel's functions. We consider here an 
effect of non-zero bandgap upon the resonance states. It is convenient to re-write Eq. ([^5)1 using the Bessel function 
explicitly (see eq (j26|) ) and choosing a — 1. We analyze here resonances for the parameters values lying in the region 
I as an example. 

Introducing the hyperbolic coordinates according to Eq. and using the dimensionless parametrization Eq. ([T7|) 
, we can write 



Jj-l/2 (0 (0 + J, + l/2 (^) i?)lV2 (^). 

- tanpi \/(e-l)/(e + l) exp {-^j^) Jj+1/2 (cr) -ffj+i/2 (^) 



r(i) 



+ v/(e + l)/(e-l)exp (V-i) J,_i/2 (0 i/)_\/2 (0 
ReE > 0, j = ±1/2, ±3/2,..., 



where — 7V — 1. If ipi — > ±cxd (asymptotical approach to the bisectrix), we have for the roots 

(0 =0, ' ±00, 



(83) 



(84) 



since the roots of the equation Hj2i/2 ("») ~ ^ ^^'^ extraneous ones. Solving Eq. (|84p . we obtain the Bessel function 
roots A^^'Fi/^)^ Therefore, 



±00 



(85) 
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This formula can be written in the form: 

= ytT^A^^^V^ (86) 

We see that appearing of non-zero mass shifts "arc's" bearing points relative to the positions, determined by Eqs.. 
([M)) . in the direction of the energy increase (see Fig. 14). This shift is essential only for small n since A,^^^^^^-* 
increases fast with n increasing. When m ^ 0, a distribution of signs of Eq. ([83]) right-hand-side for various m values 
at fixed is complicated. Arcs, in result, can be formed both clock-wise and anti-clock depending on the m value. 

We have analyzed the resonances for the region I; the dependence of the resonances distribution on the parameter 
7 = rom in other regions can be considered similarly. 

We studied in the previous section a motion of resonances along the figure-of-eight curves in the case of m = 0, 
when the parameters a, b are varied along the circles Eq. (j46p with various p values. In the non-zero gap case, we have 
found that a large enough mass destructs the figure-of-eight structure for p fixed; these figures associate into a whole 
complicated aggregate (see Fig. 15). 



V. CONCLUSION 



In conclusion, we considered the bound and resonance electronic states for the two-dimensional Dirac equation 
with the short-range perturbation. The short-range perturbation is approximated by the delta function S (r — tq) 
with different amplitudes in the upper and lower electronic bands. We have in result local perturbations of the 
potential-like and mass-like types. respectively with amplitudes a and b. 

We derived the characteristic equations for the bound and resonance states. The characteristic equation is presented 
in forms convenient both for an analysis of the bound and resonance states. Energy levels behavior in dependence on 
the perturbation amplitudes was investigated both analytically and numerically. A general picture of the electronic 
spectrum dependence on the potential, and mass local perturbation is considered for various mass and angular 
momentum values; the results are presented as diagrams in the (a, 6)-plane. Absence of a threshold for forming of the 
bound state in the vicinity of the upper bandgap edge is obtained for the lower angular momenta \j\ = 1/2 which is 
in accord with the general principles of quantum mechanics for the dimension 2-1-1. 

A qualitative distinction of > b^ and > b^ cases is shown. A countable set of eigenvalues existence domains 
in the (a, 6)-plane is present in the former case, while monotonic approaching to the asymptote takes place, when 
b^ — — > oo. in the latter one. Positions of asymptotes are determined by the parameter tanh^^ (b/a) ■ Higher 
angular momenta are investigated as well. Twisted eigenvalue domains with crossing boundaries are found. Domains 
of eigenvalues are narrowing with the angular momentum increase and at the mass value fixed. A zero-energy solution 
with vanishing imaginary part exists in the zero-gap case. This solution has to be considered as a limit at m — > 0. 

Behavior of resonances depends essentially on the sign of (a^ — 6^) . When > b^, resonance trajectories in the 
energy complex plane have a form of closed curves circumscribed periodically at moving from the coordinate origin 
to infinity in the (a, 6)-plane. Motion along the hyperbolae \/ o? — h"^ — const maps onto the countable set of arc-like 
trajectories in the energy complex plane. These arcs bear on the real axis; positions of the bearing points are exactly 
determined by the Bessel function roots. When b^ ~ a? > 0, an increase of b^ — a? forces the arcs to get flatter so that 
if — a? — > cx), the arcs asymptotically nestle up to the real axis. The resonance width is exponentially small in 
this case. For higher angular momenta j, the rate of this decrease is getting small due to the factor . Therefore, 
higher spherical harmonics do not contribute into the resonance scattering. 

Two kinds of the non-zero mass effect on the resonances behavior were found. Firstly, positions of arc-bearing 
points shift. This effect is particularly essential for lower roots. Secondly, eight-like figures resulted from mapping of 
the circular motion in the (a, 6)-plane onto the energy plane transform into large aggregates. 

The carried out analysis allowed us to describe both qualitatively and quantitatively a distribution of the bound 
and resonance states in the energy complex plane for any chosen distribution of the perturbation amplitudes a and b 
and for various magnitudes of the mass and momentum. The obtained results can be useful for understanding of the 
graphene electronic properties. 
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Figure captions 




FIG. 1: Regions of the perturbation amplitudes (a, 6)-plane. Bisectrices of tlie coordinate angles separating the (a, 6)-plane 
are presented by the thick solid lines. The singular hyperbolae 3 = 7r/2 + nn are presented by the thick dashed lines. 
Non-singular hyperbolae are presented by the thin solid lines. The hyperbolae 3 = nn are presented by the thin dotted lines. 



(E) 




FIG. 2: Schematic diagram of the arc-like resonance trajectory. En^^^^'^^ are the roots of the Bessel functions Jj±i/2 (roE) . 
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FIG. 3: Recommuting of the arc system in the region I, when crosses the singular hyperbolae p-^ = 7r/2 + nn, n = 0, 1 . . . 

(a) 7rn < < 7r/2 + Trn. 

(b) 7r/2 + vrn < Pi < tt (n + 1) . 




FIG. 4: NumericaUy calculated resonance trajectories for various p^ values in the region I. Solid lines: p^ — 0.8, dashed lines: 

Pi = 1.4, dotted line pj — 2.0. 

In the inset: Hyperbola paths in the (a, fe)-plane. 

In the body of figure: Maps of the above paths onto the resonances trajectories in the energy complex plane. 




FIG. 5: Numerically calculated resonance trajectories for various p values in the region I. Solid lines: p^ — 2.5, dashed lines: 

Pi — 3.0, dotted grey lines p^ = 3.5. 

In the inset: Hyperbola paths in the (a, 6)-plane. 

In the body of figure: Maps of the above paths onto the resonances trajectories in the energy complex plane. 
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! n ^ 2;r ReE 

2 2 

FIG. 6: Numerically calculated resonance trajectories for motion along the rays in the (a, fe)-plane in the region I. Solid lines: 
</> = 7r/7, dashed lines: = 7r/16, dotted lines: = — tt/IO. 
In the inset: Rays in the (a, b)-plane. 

In the body of figure: Resonance trajectories in the energy complex plane. 




FIG. 7: Numerically calculated resonance trajectories for various values in the region II. Solid lines: — 0-8i dashed lines: 

P2 = 1.4, dotted lines = 2.0. 

In the inset: Hyperbola paths in the (a, b)-plane. 

In the body of figure: Maps of the above paths onto the resonances trajectories in the energy complex plane. 
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FIG. 8: Numerically calculated resonance trajectories for motion along the rays in the (a, 6)-plane in the region II. Solid lines: 
4> = 57r/14, dashed lines: (p = 77r/16, dotted lines <j} = Stt/S. 
In the inset: Rays in the (a, 6)-plane. 

In the body of figure: Corresponding resonance trajectories in the energy complex plane. 
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FIG. 9: Numerically calculated resonance trajectories for motion along the circles in the (a, fo)-plane. Solid lines: p — 0.8, 
dashed lines p — 1.4, dotted lines p = 2.0. 
In the inset: Circles in the (a, 6)-plane. 

In the body of figure: Corresponding figure-of-eight resonance trajectories in the energy complex plane. 





FIG. 10: Eigenvalue dependence on the radial coordinate pj in the region I for two angle magnitudes: solid line (f) = 7r/6, 
dashed line 4> = — 7r/6.1n the inset: Two rays in the (a,fe)-plane with the above angle values. 




FIG. 11: Eigenvalue dependence on the radial coordinate p2 in the region II for three angle values (corresponding rays are 
shown in the inset): solid lines (j) — 2n/3, dashed lines (f) = it/2, dotted lines (j) — —2-k/3. 



FIG. 12: Diagram of bound states for j — 1/2 and 7 = 2. "Fronts" e = +1 are shown by solid lines, "fronts" e = —1 are shown 
by dashed lines. Thin solid line hyperbolae are asymptotes of the "fronts". 




FIG. 13: Diagram of bound states for j = 3/2 and 7 = 2. "Fronts" e = +1 are shown by solid lines, "fronts" e = —1 are shown 
by dashed lines. Thin solid line hyperbolae are asymptotes of the "fronts".. 
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FIG. 14: Arc-like resonance trajectory dependence on the mass (gap) value. Solid lines: m = 0, dashed lines: m — 0.5, dotted 
line m = 1.0. 




FIG. 15: Figure-of-eight resonance trajectory dependence on the mass (gap) value. Solid lines: m — 0, dashed lines: m — 0.25, 
dotted line m = 0.5. 



